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I. INTRODUCTION 

Precision studies, both theoretical and experimental, 
of non-relativistic QED bound-states have historically 
provided a wealth of information about fundamental 
physical parameters such as the fine structure constant 
and the masses of the electron, muon and proton [||. 
In order to keep pace with the ever-improving precision 
of current experiments, increasingly intricate theoretical 
calculations must be performed. Fortunately, the emer¬ 
gence of new calculational techniques has brought many 
such calculations within reach. 

In a previous letter |jj, we introduced a practical al¬ 
gorithm which allows a calculation of the bound-state 
energy levels in a given order of perturbation theory (in 
a and Za ) as an expansion in powers and logarithms of 
m/M with an arbitrary precision. The opposite situa¬ 
tion, the calculation of the energy levels to all orders in 
a but in a fixed order in the ratio m/M, has been studied 
in the literature H H ||, ||. We illustrated the method 
by calculating the 0(a(Za) 5 ) radiative recoil corrections 
to the average energy shift and hyperfine splitting of a 
non-relativistic QED bound-state. This analytic result 
enabled us to resolve a discrepancy between two previ¬ 
ous calculations of O(a(Za) 5 m 2 /M) corrections to the 
average energy shift, thereby removing a major source 
of theoretical uncertainty in the isotope shift (i.e. the 
difference between the 2 S to IS 1 transition energies in 
deuterium and hydrogen). 

In this paper we extend our previous results in several 
ways. For instance, we illustrate how the energy shifts 
can also be expanded as a series in (1 — m/M). Beyond its 


applicability to situations where to ~ M, this expansion 
method provides a useful cross-check on the comparably 
more difficult method of expanding in m/M. We also ex¬ 
amine the convergence properties of these expansions in 
order to ascertain some general guidelines about the ac¬ 
curacy of the truncated series that necessarily arise when 
using these expansions. In addition, we present a calcu¬ 
lation of the pure recoil corrections using an expansion 
in m/M. 

The rest of this paper is organized as follows. In 
Section II we discuss the framework of the calculation 
and our expansion method in general terms. Section III 
contains a detailed technical description of the calcula¬ 
tion. In Section IV we present and discuss the results of 
our calculations for the 0(a(Za) 5 ) radiative recoil and 
0((Za) 6 ) pure recoil corrections to the average energy 
shift and hyperfine splitting of a generic QED bound- 
state. Our results are summarized in Section V. Finally, 
the Appendix illustrates the techniques with which a 
class of loop integrals — the so-called eikonal integrals 
— can be evaluated. 


II. FRAMEWORK OF THE CALCULATION 

Non-relativistic bound-states can be conveniently de¬ 
scribed by an effective field theory which exploits a sep¬ 
aration in the energy scales to, mv, and mv 2 j7|. This 
approach is facilitated further by the use of dimensional 
regularization §B00- The basic idea is that the 
Coulomb potential supplies the dominant interaction and 
all other interactions provide corrections which can be 
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evaluated using the familiar time-independent perturba¬ 
tion theory of quantum mechanics. These corrections can 
be divided into two classes. 

The first class is that of the so-called “soft” contri¬ 
butions, governed by long-range potential terms. These 
contributions can be evaluated for arbitrary masses of the 
constituent particles because the essential soft dynam¬ 
ics of a non-relativistic bound-state, as described by the 
Schrddinger equation, are characterized by the reduced 
mass of the system rather than the individual masses of 
the consituents. As a result, once the soft contributions 
are obtained in the equal mass case 0 , the more 
general mass case follows easily. 

The “hard” contributions make up the second class, 
and these contributions lead to the interactions that can 
be characterized by S(r ) potential terms. Such terms 
result from the relativistic region of loop-momentum in¬ 
tegrals, and they are usually obtained as Taylor expan¬ 
sions of scattering amplitudes in terms of the spatial mo¬ 
mentum components of the external particles, which are 
taken to be on-shell. At lowest order in a, the hard dia¬ 
grams should be evaluated exactly at threshold, whereby 
the constituents have zero relative velocity. This implies 
that the relevant loop-momentum integrals depend on 
only two scales, m and M. The hard contributions have 
a much more complicated dependence on the mass scales 
than the soft contributions do, and this is why we will 
expand the hard scattering diagrams in powers of either 
m/M or (1 — m/M). By expanding the integrands , we are 
left with only homogeneous, one-scale integrals to evalu¬ 
ate, and this consitutes a substantial simplification. This 
method lends itself to automation so that many terms 
of the expansion can be obtained, with the only limita¬ 
tion being set by the available computing power. High- 
performance symbolic algebra software is of great help in 
such computations (we use FORM [jl2f). 

Our m/M expansion method is motivated by a proce¬ 
dure in which Feynman diagrams are expanded in large 
masses and momenta [|Ti| [l4|, |u|. Although this proce¬ 
dure was originally expressed in a different way, it can 
be reformulated more practically using the notion of mo¬ 
mentum regions. The algorithm, which is applied directly 
to the loop integrals, consists of five steps ]I(|. First, 
identify the large and small external scales in the inte¬ 
grals. Second, divide the integration volume into regions 
so that the momentum flow through any of the internal 
lines is of the order of one of the external scales. More 
specifically, the statement k ~ M asserts that k 2 > to 2 
in Euclidean space. Third, perform Taylor expansions 
within every region for any individual denominator fac¬ 
tors (propagators) where the terms within these factors 
depend differently on the external scales. Fourth, inte¬ 
grate the expanded integrands from every region over the 
initial integration volume — in other words, ignoring the 
constraints that identify the regions. Finally, add the 
contributions arising from the individual regions in or¬ 
der to obtain the final result. The fourth step of this 
algorithm requires further explanation, since by ignoring 


the constraints on the individual regions, it may appear 
as if contributions to the total integral are counted more 
than once. This does not happen because the extra con¬ 
tributions to the total integral that are introduced by 
removing the constraints on individual regions can be 
expressed as scale-less integrals, and scale-less integrals 
vanish in dimensional regularization. This implies that 
the integrals from the various regions are different ana¬ 
lytic functions of the parameters of the problem. In the 
next section, we will demonstrate this algorithm in detail. 


III. PROCEDURE 

In the previous letter [iij, we outlined the expansion 
procedure for the radiative recoil diagrams. These dia¬ 
grams are shown in Fig. [|. 
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FIG. 1: The forward-scattering radiative-recoil diagrams. 
The bold line represents the heavy constituent of the bound- 
state (e.g. a proton if we consider hydrogen) and the thin 
line — the light one (an electron). Diagrams with the crossed 
photons in the f-channel are not displayed. 

We will now describe the expansion procedure using 
the pure recoil diagrams, since a few additional compli¬ 
cations arise. These diagrams are shown in Fig. [^. 

To illustrate the method we focus on the last diagram 
in Fig. |] and consider the following scalar integral: 

/'__m 

J {kl)(kl){ki + k 2 ) 2 [(fci + fc 2 ) 2 + 2(fci + k 2 )pi] 1 

1 

X _ 

(fc 2 + 2k 1 pi)[k\ — 2k lP2 + *<5)(fcf + 2k 2 p 2 + id) 

Here [d D fc] stands for d D k/(2-7r) D , pi = rnQ , p 2 = MQ , 
where Q = (1,0, 0,0) is the time-like unit vector. Only 
the relevant infinitesimal imaginary parts of the propaga¬ 
tors have been displayed. We are going to illustrate the 
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FIG. 2: The forward-scattering pure recoil diagrams. 


expansion of the integral in Eq. (0 in powers of m/M 
following the five steps outlined above. 

There are five momentum regions to be considered. 
In the first one all the momenta are of the order of the 
large mass M. In this case one can expand the electron 
propagators in mQki. The resulting integrals are all of 
the form: 

_ [d D k 1 }[d D k 2 ] _ 

(fc2) ai (fc|) a2 (fci +fc 2 ) 2a3 (fcf - 2fc 1 p 2 )“ 4 (fc| + 2k 2 p 2 ) a s ’ 

, (2) 

with some integer powers a^. One immediately recognizes 
that all these integrals are identical with the general two- 
loop self-energy integrals of the particle with mass M for 
which the general solution is known fl7| . 

Next, there is a momentum region where fci ~ M and 
k 2 ~ rn. By using k\ ~ M to expand the electron prop¬ 
agators we obtain Eq. (j2|) again. With k 2 ~ to, the 
(fci + k 2 ) 2 and (fcf + 2 k 2 p 2 ) factors can also be expanded, 
so that the only denominator factors which depend on k 2 
are (fcf) and (2 k 2 p 2 ). Since we are working at threshold, 
we have p 2 = MQ , which leads to 

/'_M_ = 0 (3) 

J (k 2 ) a (2k 2 p 2 + i6)P ’ l J 

so that this region provides no contribution to the am¬ 
plitude for this particular diagram. 

The third momentum region has k\ ~ M and k 2 ~ M 
but k\ + k 2 ~ to. After a Taylor expansion in small vari¬ 
ables, the integrals in this region factorize into products 
of two simple one-loop integrals. 


In the fourth region, ki ~ to and k 2 ~ M. A Taylor 
expansion in small variables allows the integrals in this 
region to be factored into one-loop integrals as 


[d^r] 


(. kf) ai (kf + 2kipi)° 2 (2k\pi) a3 

f ^ [d c fc 2 ] 

J (fc|) a4 (fc| + 2/c 2 p 2 )“ 5 ’ 


(4) 


The fc 2 integral is a trivial one-loop integral. The k\ 
integral can be converted to the same simple form, along 
with integrals like Eq. <® , by multiplying it by factors of 

(kl + 2kip 1 ) (2k\p\) 

(A?) (fc?) 1 j 


until either a 2 or a 3 is brought to zero. 

The fifth region is characterized by the condition ~ 
k 2 ~ to. In this case, the heavy particle propagators can 
be expanded into static, or as we will call them, eikonal, 
propagators. The integrals in this region are of the form 


[d D fci][d°fc 2 ] 


{kl) ai {kl) a2 (k\ + fc 2 ) 2 “ 3 [(fei + k 2 ) 2 + 2(ki + k 2 )pi] 


(fc 2 + 2 fcipi) a 5 ( 2 fcip 2 — iS) ae (2k 2 p 2 + iS ) a7 


( 6 ) 


Notice how the eikonal propagators arising from ( k 2 — 
2 kp + iS) factors acquire — i6 pole terms. Such terms are 
important in this region and must be carefully accounted 
for. The integral in Eq. © can be simplified using the 
identity 


1 1/1 1 

(fc 2 + 2kp)[2kp) k 2 \( 2 fcp) (k 2 + 2kp) 

Once one of the seven factors in Eq. © has been re¬ 
moved, an identity can be constructed from the observa¬ 
tion that the remaining six factors are linearly dependent. 
Using such an identity, the integrals in this fifth region 
can be expressed as one of four types of integrals. The 
first type is the two-loop self energy integral of a particle 
with mass to for which the general solution is known 0 - 
The remaining three types contain eikonal propagators: 



Ef = 


E? = 


Et* = 


[d D fci][d°fc 2 ] 


(fc 2 ) ai (fc |)“ 2 (fci — k 2 ) 2a3 
1 

+ 2k 2 pi) ai (2kipi ± z<5) a s ’ 

f [d^fcrHd 2 ^] 

(fc 2 ) ai )“ 2 (k 2 + 2kipi) 2as 

1 

[(fei + k 2 ) 2 + 2(ki + k 2 )pi] ai (2k 2 pi ± iS) a5 
f [d D /d][d D k 2 ] 


( 8 ) 


(9) 


{ki) ai (Jil) 02 {k\ — k 2 ) 2a3 

1 

(2k\p 2 ± iS) a4 (2k 2 p 2 ± iS ) a5 


( 10 ) 
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At threshold, the £3 integrals are exactly zero for similar 
reasons as are needed to establish Eq. (||). This leaves 
us with two new types of integrals, £1 and £ 2 , required 
for the recoil calculations. The calculations for radiative 
recoil diagrams involve £2 but not £1 integrals. In Ap¬ 
pendix A, we outline the procedure by which the £1 and 
£2 integrals can be evaluated. 

We shall now describe a second method of expansion, 
relevant to the scenario where the two bound-state con¬ 
stituents have similar, but not necessarily equal, masses. 
Although this scenario is not realized by any common 
QED bound-states, this second expansion can provide a 
useful check on the first expansion method. 

The essential idea is to introduce an expansion param¬ 
eter 


only diagrams we have to consider. We have done the 
calculation in a general covariant gauge; the cancellation 
of the gauge parameter dependence serves as a check of 
the computation. 

For the S -wave ground state energy £ we define 

£ = £aver + f ^ - SjQ j £hfs , (14) 


where J = 0,1 is the total spin of the two fermions form¬ 
ing the bound-state. 

For the hyperfine splitting we obtain 




8 (Za) 4 p 3 
3 mM 


a(Za) 



TO 


( 11 ) 


so that the external momentum of the light particle, pi, 
can be written in terms of y and the external momentum 
of the heavy particle, p 2 , via 


Pi = (1 - y)P2 • (12) 


Then, any massive propagator containing p i can be ex¬ 
panded, as a series in powers of y, in terms of the corre¬ 
sponding propagators containing p 2 : 


1 

k 2 + 2 kpi 


v+ (2ykp 2 r 

^ ( k 2 + 2 kp 2 ) n + 1 ' 


(13) 


As a result, the two-scale general scalar integral in Eq. ( 1 ) 
is expanded, as a series in powers of y, in terms of two- 
loop on-shell self-energy integrals. As an additional cross¬ 
check, we note that in the limit that the two masses are 
equal, we have y = 0 , so that the positronium results 
|j], [ToJ can be recovered from the leading term of this 
expansion. 

Another method has been proposed [jl 8 | for deal¬ 
ing with similar problems involving more than one 
mass/energy/distance scale. Intermediate parameters 
are introduced to separate various scales and the calcu¬ 
lations are performed in four dimensions. An advantage 
of such an approach is that it avoids various complica¬ 
tions arising when working in D dimensions. However, 
it spoils the homogeneity of integrals and it is not clear 
whether one can apply integration-by-parts algorithms, 
which are crucial for larger calculations involving many 
terms of expansions. 


IV. RESULTS AND DISCUSSION 

We have applied our algorithms to compute the 
0(a(Za) 5 ) radiative recoil corrections to the average en¬ 
ergy shift and the hyperfine splitting of a general QED 
bound-state composed of two spin- 1/2 particles with 
masses to and M. In this case the soft contribution is 
absent and the hard corrections shown in Fig. 1 are the 
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where p = mM /(m + M) is the reduced mass of the 
bound-state. 

For the spin-independent energy shift we find 
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To our knowledge the terms 0(rn 3 /M 3 ) and higher 
are new for both E^{ s and E avf , r , while the other terms 
have been obtained previously |lfj] . The coefficient of the 
0{m/M) term in Eq. @ was the subject of some con¬ 
troversy, since two different numerical results have been 
reported, [^ 0 ], ^ 1 ], ^ 2 | and J23|. 

Our result for this term, 


y(Za) 


5 P 


i 2 M 


^-^-21n2 


~ -1.32402796 a(Za) 


5 M 3 m 

m 2 M ’ 


(17) 


is in excellent agreement with the numerical result of 
Ref. Q where the coefficient —1.324029(2) was ob¬ 
tained, and has since been confirmed in an independent 
analytical calculation |^| . 

The expansions in m/M of Eqs. ( |l5| ) and (jl^) do not 
yield accurate results for increasing values of m/M , even 
though an untruncated series can be expected to con¬ 
verge on the interval m/M £ [0,1). The convergence of 
the terms in the series which we have calculated is de¬ 
picted in Fig. 3. The upper graph shows the hyperfine 
splitting calculations as a function of m/M and the lower 
graph shows the corresponding results for the average en¬ 
ergy shift. In both cases, several curves are plotted, each 
representing the sum of the first N terms of the expan¬ 
sion, where N is shown in the legend. The graphs suggest 
that more terms in the series would be required to obtain 
reliable values of the hyperfine splitting for m/M larger 
than about 0 . 2 ; our series for the spin-independent en¬ 
ergy shift should be reliable for m/M values up to about 
0.5. 

The expansion in y = 1 — m/M aims to address the 
region where m/M is large. For the hyperfine splitting 
we find 
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FIG. 3: 0(a(Za ) 5 ) radiative recoil contributions to the hy¬ 
perfine splitting and average energy shift in the m/M expan¬ 
sion. 
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In Fig. [| we illustrate the convergence of these series. 
Both graphs suggest that the terms we have calculated 
in this expansion should yield reliable results for m/M 
larger than about 0.15. 
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FIG. 5: 0(a(Za) 5 ) radiative recoil contributions to the liy- 
perfine splitting and average energy shift for arbitrary values 
of m and M. 
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FIG. 4: 0(a(Za) 5 ) radiative recoil contributions to the hy- 
perfine splitting and average energy shift in the 1 — m/M 
expansion. 

Although neither expansion can handle arbitrary val¬ 
ues of to and M, the bound-state energy level correc¬ 
tions for any m/M ratio can be reliably calculated with 
one of the expansions. To illustrate this, we have spliced 
together the two expansions in Fig. 5. The expansions 
merge nicely at an m/M value of 0.15, thereby providing 
a useful check on these methods. It is also important, in 
the context of future applications of these methods, to 
have covered the entire range of m/M. 

We have also calculated the corresponding 0((Za) 6 ) 
pure recoil corrections. For the hyperfine splitting we 
obtain 
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The divergences in this result are canceled by soft-scale 
terms, which can be calculated by extending the calcula¬ 
tion of 10 to the unequal mass case, resulting in 
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Combining the results of Eqs. (|2^) and (^l|), we find that 
the total (Zay pure recoil contribution to the hyperfine 
splitting is 
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The terms in the first line of Eq. (22) are in agreement 
with the result first obtained in [25]. The remaining 
terms, arising solely from the hard-scale contributions in 
Eq. (20), can be used to obtain an analytic approxima¬ 
tion to the function /( x) near x = 1 in Eq. (72) of |f2tS| . 

For the hard contribution to the spin-independent en¬ 
ergy shift we find 
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The m/M term of this expansion is in agreement with a 
calculation in Ref. Q. To our knowledge, the subsequent 
terms of this expansion are new. In addition, we have 
calculated these energy level shifts as an expansion in 
(1 — m/M), but for brevity we shall omit these results. 

In spite of the fact that the hard-scale contribution 
to the average energy shift given by Eq. (^3|) is finite 
at this order, soft contributions are also present and are 
needed to arrive at the physical result for this quantity. 
These soft contributions can be obtained by a calculation 
completely analagous to the one that produced Eq. (21). 


V. CONCLUSIONS 


it might be feasible to extract the terms enhanced by 
one or more factors of ln(A 1/m) by examining the sin¬ 
gularities of the contributions from different expansion 
regions. Since these singularities must cancel in the com¬ 
plete result, their coefficients can be found by a partial 
calculation of the divergent parts of those contributions 
which can be evaluated most easily. 

In a more general context, these expansion techniques 
are applicable to a plethora of other types of problems in¬ 
volving multiloop calculations with more than one exter¬ 
nal scale. Many kinds of superficially disparate physical 
problems often depend on a few common classes of loop 
integrals, thereby reducing the number of technical hur¬ 
dles which restrict the progress of precision calculations 
in particle physics. 
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Appendix A: Calculation of Eikonal Integrals 


An expansion in powers of m/M typically gives 
rise to integrals containing eikonal propagators such as 
(2 kpzkiS), arising from expansions in the momentum re¬ 
gion where the loop momenta are all ~ to. The easiest 
way to solve them is to employ the integration-by-parts 
techniques [ fl7[ |28[| so that any integral of the form in 
Eqs. (jS]) and(P) can be algebraically expressed as a com¬ 
bination of the two-loop on-shell self-energy integrals and 
four new master integrals. The latter are the only inte¬ 
grals we have to compute, and the results read 
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(An) 0 


[d D k 1 ][d D k 2 


I ( k\Q — 1 ± *<5)(fc| + iS) [(A* + k 2 ) 2 — 1 + id] 
[2T(1 - e)T(3e - 2)R(4e - 3, 2e - 1) 


-(l T l)v^r { 2e — 




(24) 


We have demonstrated a method by which the correc¬ 
tions to the energy levels of a QED bound-state, with 
constituents of mass to and M , can be expanded in ei¬ 
ther powers of m/M or (1 — m/M). Both expansions are 
applied directly to the integrands of the loop integrals 
arising from the hard-scale contributions to the energy 
shifts. We have demonstrated the utility of these proce¬ 
dures by computing several terms in the expansions for 
the a(Za) 5 radiative recoil and ( Za) e pure recoil correc¬ 
tions to both the average energy shift and the hyperfine 
splitting of a general QED bound-state. 

Further studies of QED bound-state problems, us¬ 
ing the methods described in this paper, might involve 
higher-order corrections to the energy level shifts. Even 
in the absence of a complete calculation of such terms, 


4 = 


[d D fci][d D fc 2 


= ± 


(kiQ ± iS)(k 2 — 1 + iS) [(fci + k 2 ) 2 — 1 + iS\ 
3 


(4tt) 


D 


T 2e- 


B \ -- + £,—- + £ )(25) 


Please note that there is a typographical error for 
in Q. For clarity, we now outline the process by which 
these master integrals are calculated. 

Starting with j //, we can Wick rotate the momenta 
and subsequently ignore the poles. After shifting k\ to 
k\ — k 2 and using the identity 


1 1 f A ^ -1 

~ B(a, (3) J 0 dX [A + B\] a +0 ’ 


( 26 ) 




















we have 

J+ = 


[i d D ki][d D k2] 


I 


k$ 


dX 


o [^1 + 1 + ^{kiQ — feQ + 1 )]^ 


■ (27) 


The ki integral can be evaluated, after completing the 
square and using Q 2 = — 1 , so that 


J+ = 


r(e) 

(4 tt) d / 2 


dX 


[d D k 


k 2 —b A + 1 — Xk 2 Q) 


. (28) 

Introducing a second parameter from the identity in 
Eq. (p(i|), the k 2 integral can be evaluated, yielding 


Jt = 


T(2e — 1) 
(47r) 


D 


x dp p 

J o 

With the substitution 


dX 

-e f P^_ 

4 


+ A + l 


1 — 2 e 


p = 4z 


(A/2 + 1 ) 2 
A 2 


(29) 


(30) 


we can readily integrate over the remaining parameters 
to obtain the result for J+ in Eq. ( p4| ) . The evaluation of 
j/~ starts by introducing a Feynman parameter to com¬ 
bine and integrate over the A^-dependent factors, and 
is followed by the introduction of a parameter from the 
identity in Eq. (26). 

Turning to Jj - , we note that the —id pole term in the 
eikonal propagator prohibits a Wick rotation. Instead, 
the identity 


dx 


= P 


dx 


Loo x-a 


+ inS(x — a) , (31) 


J-oo x — a±iS 
with respect to the h® integral, allows us to write 

J- = J+ + AJi (32) 


with 


AJi = 27ri 


[d D k\]\d D k 2 ] S(k° — 1) 

(fc| + id)[(Aq + k 2 ) 2 — 1 + id] 


(33) 


A Feynman parameter can be introduced to evaluate the 
k 2 integral in AJi. After performing the A;J integral, we 
have 


AJi = 


-27rr(e) 


dx 


(47r) I3 / 2 (27r) 13 jQ 

f d D ~ 1 k 


[x — *(1 — x) + x{l — a;)k 2 ] £ 


, (34) 


where k denotes the (D — l)-dimensional spatial momen¬ 
tum associated with ki. Working in hyperspherical coor¬ 
dinates, and using 


iD-l 


n = 


we obtain 


2 tt ( d - 1 )/ 2 

WT ’ 


(35) 


AJi = - 


^ r(e) f 
(4ir) D r (| - e) J 0 
r°° k D ~ 2 dk 


dx 


I o [x 2 + x(l — x)k 2 ] c 


with k = Ikl. With the substitution 


(36) 


k = 


1 — x 


(37) 


the x and z integrals can be evaluated in terms of B- 
functions, yielding the result in Eq. (pl[). The evaluation 
of J 2 proceeds in a similar fashion. 
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